/************************************************************************
 * MechSys - Open Library for Mechanical Systems                        *
 * Copyright (C) 2012 Sergio Torres, Diego Fregolente                   *
 *                                                                      *
 * This program is free software: you can redistribute it and/or modify *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * any later version.                                                   *
 *                                                                      *
 * This program is distributed in the hope that it will be useful,      *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the         *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with this program. If not, see <http://www.gnu.org/licenses/>  *
 ************************************************************************/


// Std lib
#include <math.h>

// MechSys
#include <mechsys/dem/domain.h>
#include <mechsys/util/fatal.h>
#include <mechsys/mesh/unstructured.h>
#include <mechsys/mesh/structured.h>

using std::cout;
using std::endl;
using std::ofstream;
using DEM::Domain;
using Util::PI;


int main(int argc, char **argv) try
{
    
    DEM::Domain dom;

    if (argc<2) throw new Fatal("This program must be called with one argument: the name of the data input file without the '.inp' suffix.\nExample:\t %s filekey\n",argv[0]);
    
    //number of threads
    size_t Nproc = 4; 
    if (argc==3) Nproc=atoi(argv[2]);
    String filekey  (argv[1]);
    String filename (filekey+".inp");
    if (!Util::FileExists(filename)) throw new Fatal("File <%s> not found",filename.CStr());
    ifstream infile(filename.CStr());
        

     //set the simulation domain ////////////////////////////////////////////////////////////////////////////
    //
    
    
    double verlet;      // Verlet distance for optimization
    String ptype;       // Particle type 
    size_t VOut;        // Decide is video should be render
    double fraction;    // Fraction of particles to be generated
    double Kn;          // Normal stiffness
    double Kt;          // Tangential stiffness
    double Gn;          // Normal dissipative coefficient
    double Gt;          // Tangential dissipative coefficient
    double Mu;          // Microscopic friction coefficient
    double R;           // Spheroradius
    size_t seed;        // Seed of the ramdon generator
    double dt;          // Time step
    double dtOut;       // Time step for output
    double Lx;          // Lx
    double Ly;          // Ly
    double Lz;          // Lz
    size_t nx;          // nx
    size_t ny;          // ny
    size_t nz;          // nz
    double rho;         // rho
    double Tf;          // Final time for the test
    {
        infile >> verlet;       infile.ignore(200,'\n');
        infile >> ptype;        infile.ignore(200,'\n');
        infile >> VOut;         infile.ignore(200,'\n');
        infile >> fraction;     infile.ignore(200,'\n');
        infile >> Kn;           infile.ignore(200,'\n');
        infile >> Kt;           infile.ignore(200,'\n');
        infile >> Gn;           infile.ignore(200,'\n');
        infile >> Gt;           infile.ignore(200,'\n');
        infile >> Mu;           infile.ignore(200,'\n');
        infile >> R;            infile.ignore(200,'\n');
        infile >> seed;         infile.ignore(200,'\n');
        infile >> dt;           infile.ignore(200,'\n');
        infile >> dtOut;        infile.ignore(200,'\n');
        infile >> Lx;           infile.ignore(200,'\n');
        infile >> Ly;           infile.ignore(200,'\n');
        infile >> Lz;           infile.ignore(200,'\n');
        infile >> nx;           infile.ignore(200,'\n');
        infile >> ny;           infile.ignore(200,'\n');
        infile >> nz;           infile.ignore(200,'\n');
        infile >> rho;          infile.ignore(200,'\n');
        infile >> Tf;           infile.ignore(200,'\n');
    }


//    Domain d;
//    d.CamPos = Vec3_t(0.0, 0.0, 90.0); // position of camera
    double cam_x=0.0, cam_y=340.0, cam_z=0.0;
    dom.CamPos = cam_x, cam_y, cam_z;


    // particle
    bool Cohesion = false;
    if      (ptype=="sphere")  dom.GenSpheres  (-1, Lx, nx, rho, "HCP", seed, fraction);
    else if (ptype=="voronoi") dom.AddVoroPack (-1, R, Lx,Ly,Lz, nx,ny,nz, rho, Cohesion, !Cohesion, seed, fraction);
    else if (ptype=="tetra")
    {
        Mesh::Unstructured mesh(/*NDim*/3);
        mesh.GenBox  (/*O2*/false,/*V*/Lx*Ly*Lz/(0.5*nx*ny*nz),Lx,Ly,Lz);
        dom.GenFromMesh (mesh,/*R*/R,/*rho*/rho,Cohesion,false);
    }
    else if (ptype=="rice") dom.GenRice(-1,Lx,nx,R,rho,seed,fraction);
    else throw new Fatal("Packing for particle type not implemented yet");
    
    dom.Center(Vec3_t(0.0,0.0,110));
    dom.Alpha = verlet;

     

    dom.AddSphere(-4, Vec3_t (0,0,75), 3, 10*rho);                                                          ///< Add sphere
    
    
    Vec3_t X2(0,0,47);    
    Vec3_t X3(0,0,180);    

    
    dom.AddCylinder (/*int Tag*/-3, /*Vec3_t const & X0*/X2, /*double R0*/13.5, /*Vec3_t const & X1*/X3, /*double R1*/60, /*double R*/R, /*double rho*/200*rho);       

    Vec3_t X02(0,0,47);    
    Vec3_t X03(0,0,13);    
    
    
    dom.AddCylinder (/*int Tag*/-7, /*Vec3_t const & X0*/X02, /*double R0*/14, /*Vec3_t const & X1*/X03, /*double R1*/21.5, /*double R*/R, /*double rho*/200*rho);       
    
    
    Vec3_t X0(0,0,13);    
    Vec3_t X1(0,0,-80);    
    
    dom.AddCylinder (/*int Tag*/-5, /*Vec3_t const & X0*/X0, /*double R0*/22, /*Vec3_t const & X1*/X1, /*double R1*/22, /*double R*/R, /*double rho*/200*rho);  
    dom.AddPlane(/*tag*/-6,/*position*/0.95*X1,/*spheroradius*/R,/*Lx*/100,/*Ly*/100,/*rho*/200*rho);

    dom.AddRice(-8 ,Vec3_t (0.0,0.0,40.0), 1.5,10, 200*rho,M_PI/2,&OrthoSys::e0); 

    dom.AddRice(-9 ,Vec3_t (0.0,0.0,18.0), 1.5,10, 200*rho,M_PI/2,&OrthoSys::e0);      
    
    dom.AddRice(-10,Vec3_t (10.0,0.0,18.0), 1.2,10, 200*rho,M_PI/2,&OrthoSys::e0);      

    dom.AddRice(-11,Vec3_t (-10.0,0.0,18.0), 1.2,10, 200*rho,M_PI/2,&OrthoSys::e0);      

    dom.AddRice(-12,Vec3_t (0.0,0.0,29.0), 1.2,15, 200*rho,M_PI/2,&OrthoSys::e1);
     
    
    Dict B;
    B.Set(-1 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,Mu );
    B.Set(-3 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-4 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-5 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-6 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-7 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-8 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-9 ,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-10,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-11,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    B.Set(-12,"Kn Kt Gn Gt Mu",Kn,Kt,Gn,Gt,0.0);
    dom.SetProps(B);

    
    //Fix the obstacles
    
    dom.GetParticle(-3 ,true)->FixVeloc();
    dom.GetParticle(-4 ,true)->FixVeloc();
    dom.GetParticle(-5 ,true)->FixVeloc();
    dom.GetParticle(-6 ,true)->FixVeloc();
    dom.GetParticle(-7 ,true)->FixVeloc();
    dom.GetParticle(-8 ,true)->FixVeloc();
    dom.GetParticle(-9 ,true)->FixVeloc();
    dom.GetParticle(-10,true)->FixVeloc();
    dom.GetParticle(-11,true)->FixVeloc();
    dom.GetParticle(-12,true)->FixVeloc();
    
    //Initialize the gravity on the particles
    for (size_t i=0;i<dom.Particles.Size();i++)
    {
        dom.Particles[i]->Ff = dom.Particles[i]->Props.m*Vec3_t(0.0,0.0,-9.8);
    }

    dom.Solve     (/*tf*/Tf, /*dt*/dt, /*dtOut*/dtOut, NULL, NULL, filekey.CStr(),VOut,Nproc);
    dom.Save(filekey.CStr());

    return 0;
}
MECHSYS_CATCH